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Abstract 

This  paper  deals  with  the  derivation  of  equations  suitable  for  the 
computation  of  elastic  curves  on  the  sphere.  To  this  end  equations 
for  the  main  invariants  of  spherical  elastic  curves  are  given.  A  new 
method  for  solving  geometrically  constraint  differential  equations  is 
used  to  compute  the  curves  for  given  initial  values.  A  classification  of 
the  fundamental  forms  of  the  curves  is  presented. 

1      Introduction 

Two  current  trends  in  Geometric  Modeling  are  concerned  with 

•  the  development  of  spline  techniques  on  surfaces  (see  [14],  [12],  [10]) 

•  the  use  of  nonlinear  spline  curves  of  minimal  elastic  energy  for  the 
modeling  of  smooth  shapes  (see  [2],  [3],  [4]). 

This  paper  blends  both  topics  in  considering  elastic  curves  on  the  sphere. 

In  euclidean  space  an  elastic  curve  can  be  viewed  as  an  arc  length  para- 
metrized "cubic"  spline  in  tension,  i.e.  an  elastic  curve  is  a  critical  point  of 
the  functional 


<I> 


{y)  =  /    <y",y"  > +v  <y\y' >  ds  .(1) 

Jo 


in  the  space  F  of  smooth  maps 

y:[0,/]-»R",  \y'\  =  l 

y(0)  =  P0,  y(l)  =  Pi,  y'(0)  -  V&,  y'(l)  =  Vx 

where  P0,  Px  €  Rn,  V0  6  TPoRn,  V7!  €  TPlRn,  cr  €  R  are  fixed  and  /  is 
variable.  Expressing  y"  in  the  Frenet  frame  yields  the  functional  $  in  the 
form 

*(y)=   l\<s)Y  +  ads  (2) 

which  explains  the  interest  in  elastic  curves  as  those  curves  that  minimize 
bending. 

The  notion  of  cubic  splines  can  be  generalized  to  curves  on  a  Riemannian 
manifold  M  by  replacing  the  usual  derivative  of  the  tangent  vector  field 
y'  by  the  covariant  derivative  compatible  with  the  metric  of  M  (see  [13]). 
Generalizing  the  functional  (1)  in  this  way,  one  obtains  the  concept  of  elastic 
curves  on  arbitrary  manifolds.  In  the  case  of  surfaces  embedded  in  R3  the 
algebraic  value  of  the  covariant  derivative  of  the  tangent  vector  field  y'  of  a 
surface  curve  y  is  called  the  geodesic  curvature  Kg  of  y  (see  [5]  and  section 
2).  Therefore  we  may  define  an  elastic  curve  on  a  surface  S  :  A  C  R2  — ■»  R3 
as  an  extremal  point  of  the  functional 


$ 


(y)=    [\K3(s)y  +  crds,  (3) 

Jo 


in  the  subset  F  of  F  formed  by  curves  lying  on  S. 

In  section  3  a  set  of  differential  equations  for  elastic  curves  on  the  sphere  is 
derived.  This  set  includes  a  differential  equation  for  the  geodesic  curvature  of 
spherical  elastica.  Since  the  normal  curvature  of  a  spherical  curve  is  constant, 
the  differential  equation  for  the  geodesic  curvature  suffices  to  compute  the 
ordinary  curvature  of  a  spherical  elastica.  Furthermore,  a  formula  is  given 
that  expresses  the  squared  torsion  of  a  spherical  elastica  as  a  rational  function 
of  its  curvature. 

In  section  4  we  describe  the  numerical  algorithm  used  to  integrate  the 
set  of  differential  equations  derived  in  section  3.  The  equations  have  a  very 
particular  structure  defined  by  a  number  of  constants  of  motion,  and  in 
particular  they  constrain  the  elastic  curves  to  lie  on  the  sphere.  We  employ 
an  algorithm  introduced  by  Crouch  and  Grossman  [8]  which  preserves  the 
constraints  exactly. 

Since  Euler's  fundamental  work  on  plane  elastic  curves  it  is  known  that 
these  curves  can  be  classified  according  to  their  shape.  The  tools  developed 
in  this  paper  enable  us  to  present  the  fundamental  forms  of  spherical  elastica 
in  the  last  section. 


2      Geometric  Preliminaries 

Let  S  :  U  C  R  — *  R3  denote  a  regular  parametric  surface  and  N  the  unit 
normal  vector  field  of  S.  A  curve  x  :  I  C  R  — *  R3  is  a  curve  on  the  surface 
5  if  and  only  if  x  =  5  o  c  where  c  :  /  — >  (7  is  a  plane  curve  in  U.  The  unit 
normal  of  S  along  a  surface  curve  x  will  be  denoted  by  n  :=  N  o  c. 

The  Darboux  frame  &i,&2>&3  along  x  is  the  orthogonal  frame  defined  by 

bi(t)  =  |p|j|r,      62(0  =  »(*)  x  &i(t),      63(i)  =  n(t). 

The  equations  that  express  the  derivatives  61,62,63  in  the  Darboux  basis 
&!,  62,  ^3  are  given  by: 


6':       =      lL>K5&2  +^n&3,  (4) 

&2       =       -^K56!  +CJT5&3,  (5) 

63    =    -ujKnb1-uTgb2  (6) 

with  w(*)  =  |2:'(t)|. 

The  functions  Kg.  Kn  and  rg  are  called  geodesic  curvature,  normal  cur- 
vature and  geodesic  torsion.  The  geodesic  torsion  and  the  absolute  value  of 
geodesic  and  normal  curvature  are  invariant  under  reparametrization  of  the 
surface. 

The  geodesic  curvature  of  a  surface  curve  x  at  a  point  x(t)  is  the  ordinary 
curvature  of  the  plane  curve  generated  by  orthogonal  projection  of  x  onto 
the  tangent  plane  of  S  at  x(t).  It  can  be  computed  using  the  formula 


k„  = 


[x',x",  n) 


9  '         \x'P 


(7) 


A  surface  curve  with  identically  vanishing  geodesic  curvature  is  called  a 
geodesic  of  the  surface. 

The  absolute  value  of  the  normal  curvature  of  x  at  a  point  x(t)  is  the 
curvature  of  the  intersection  of  S  with  the  plane  through  x(t)  spanned  by 
the  vectors  x'(t)  and  n(t).  While  the  geodesic  curvature  is  the  curvature  of  a 
surface  curve  from  a  viewpoint  in  the  surface,  normal  curvature  measures  the 
curvature  of  the  curve  that  is  due  to  the  curvature  of  the  underlying  surface. 
If  k  denotes  the  ordinary  curvature  of  the  space  curve  x  the  identity 

k2  =  <  +  *1  (8) 

holds. 


The  geodesic  torsion  of  a  surface  curve  x  at  a  point  x(t)  is  the  torsion 
of  the  geodesic  that  meets  x  at  x(t)  with  common  tangent  direction.  A 
curvature  line  of  x,  i.e.  a  curve  with  a  tangent  vector  that  points  into  one  of 
the  principal  directions  of  the  surface,  is  characterized  by  vanishing  geodesic 
torsion. 

3      The  differential  equations  of  spherical  elas- 
tica 

Let  5  be  a  parametrization  of  a  patch  on  a  sphere  S2  C  R3  of  radius  r  and 
center  0  such  that 

S  =  rN 

and  x  be  an  arc  length  parametrized  (i.e.  |x'|  =  1)  curve  on  S.  Then 

x   =  rn   =  rb3. 
In  this  situation  (6)  implies  that 

Kn      = ,  Tg       =      0.  (9) 

r 

Since  the  absolute  value  of  Kn  and  r5  are  invariant  under  reparametrization, 
these  equations  imply  that  any  spherical  curve  is  a  curvature  line  with  con- 
stant normal  curvature. 

We  consider  the  variational  problem  of  minimizing 

/  (k5(s))2  +  a  ds 
Jo 

in  the  space  F  of  C°°  smooth  maps 

r,:[0,/]-S2,  1^1  =  1 

y(0)  =  Po,  !/(l)  =  Pu  y'{0)  =  V0,  y\\)  =  Vx 

where  P0,  Pi  £  S2,  Vo  G  Tp0S2,  V\  G  Tp^2,  cr  G  R  are  fixed  and  /  is 
variable. 

From  (4)  one  obtains  the  relation 


rcia«2  +  . 


1 

so  that  we  wish  to  minimize  the  functional 

rl\b[(s)\2  +  6ds 

■    -4 


where 

*  =  ,-! 

under  the  constraints 

IM2  =  1)      h  =  x'i      \x\ 
Hence,  we  can  apply  the  Euler- Lagrange  equations  to  the  functional 

F  =  16'J2  +  6  +  A(|61|2  -  1)  +  ii(\x\2  -  r2)  +  2  <  A,x'  -  b,  > 

to  obtain  the  differential  equations  which  govern  the  extremals: 

fix-A'  =  Q  (10) 

A&!-&;'  =  A.  (11) 

Combining  these  equations  yields 

A'fcj  +  \b[  -b'C  =  fix.  (12) 

The  derivatives  of  6}  expressed  in  the  Darboux  basis  are  given  by: 

b[     =     Kgb2 63  (13) 

K    =     *;*2  -  («?  +  ^)*i  (14) 

V"    -    (-3/c5k')6i  +  (/c'/-/Cj -— «5)62 


1,  ,       1 


+  ;K  +  ^)&3.  (15) 

Substituting  these  derivatives  and  x  =  r&3  into  (12)  and  rearranging  gives 

(A'  +  3K3K/)b1 


3       1 


+     (Ak5 -<  +  /c   +— «5)6 


r2 

2        1 


2 


-     (7  +  /t  +  ;K  +  -))63  =  0. 
Finally,  the  linear  independence  of  the  vectors  6l5  62  and  63  implies  that 

\  =  J-k)  +  C  (16) 

and 


In  order  to  determine  the  constant  C  in  terms  of  the  tension  parameter 
cr,  we  consider  the  boundary  condition 

F  <'>  -  d£(l)Al)  -  S(,)6'('> = ° 

for  the  extremal  x.  This  condition  is  implied  by  the  fact  that  the  total  length 
/  of  the  curve  is  variable  in  the  variation  (see,  e.g.,  [1]).  Thus, 

-  *J(0  -  p  + s  - 2  <  MO, *'(0  >=  o-  (is) 

Substituting  A  according  to  (11)  into  the  scalar  product  <  A,x'  >  yields 

<A(/), *'(/)>     -    <\bl{l)X{l)>-<h'{{l)Ml)> 
=    A  +  (k3(1))2  +  I 

=  -4k«))2+^+c. 

Substituting  this  expression  for  the  scalar  product  into  (18)  we  obtain 

c  +  i-^t-?)-  (19) 

These  results  are  summarized  in  the  following  theorem. 

Theorem  1  An  elastic  curve  x  under  tension  a  on  the  sphere  of  radius  r 
satifies  the  differential  equations 

x[\        [      0         1/r      0    \ 


x'n  -1/r       0       K 


3 


(20) 

0  -Kg  0       J 

where  x\  =  x,  x-i  =  rx' ,  x$  =  x  x  x'  and  where  the  geodesic  curvature  Kg  of 
x  is  a  solution  of 

<+^  +  (^-f  )«.  =  «■  (21) 

The  curvature  of  a  spherical  elastic  curve  can  be  obtained  from  (8)  and  the 
fact  that  the  normal  curvature  of  spherical  curves  is  constant.  The  squared 
torsion  of  a  spherical  elastica  can  be  expressed  as  a  rational  function  of  its 
curvature. 

Theorem  2    The  curvature  k  and  the  torsion  r  of  a  spherical  elastic  curve 
obey  the  relation: 

rW--!<«-!(± -,)««  +  I(f -§)  +  <*. 


Proof:  For  the  invariants  k  and  r  of  an  arc  length  parametrized  curve  x  in 
R3  the  relation 

r/c2  =  [&!,&;,&'/]  (22) 

holds  where  [a,b,c]  denotes  the  determinant  of  three  vectors  in  R3. 

For  a  spherical  curve  b\  and  b"  can  be  expressed  in  the  Darboux  basis  as 
follows: 

&'     =     Kgb2 b3 

r 

K  =  -(£  +  «J)*i  +  <»»■ 

Substituting  the  derivatives  into  (22)  and  squaring  yields 


1 

7' 


Since  the  differential  equation  (21)  can  be  integrated  to 

ft)1  =  ft  -  i«;  -  £  -  f  )«2, 

one  obtains  the  claimed  equation  using  (8). 

4      Tracking  elastic  curves  on  the  sphere 

The  problem  we  consider  here  is  that  of  numerically  integrating  equations 
(20)  and  (21)  in  Theorem  1.  One  can  of  course  simply  integrate  the  equa- 
tions using  a  standard  numerical  package,  such  as  an  IMSL  Runge  Kutta 
routine.  However,  the  system  of  equations  possesses  a  very  special  structure. 
As  pointed  out  in  [3],  equation  (21)  may  be  integrated  directly  in  terms  of 
Jacobi's  elliptic  functions.  We  give  more  details  of  this  process  in  the  next 
section  where  we  classify  the  various  extremals.  As  for  equation  (20)  we 
note  that  the  components  of  the  state  vector  \x\ ,x^,  xJ]T  satisfy  algebraic 
constraints  consistent  with  the  fact  that  the  matrix  [a-!,  x2i  £3]  is  simply  a 
multiple  r,  of  a  rotation  matrix.  When  a  standard  integration  package  is 
applied  to  the  set  of  differential  equations  (20)  and  (21),  these  constraints 
are  not  preserved  exactly,  and  in  particular  the  norm  of  the  vector  x^  will 
not  remain  at  the  constant  value  r.  This  is  a  particularly  important  fact 
when  we  wish  to  integrate  the  equations  over  a  large  number  of  time  steps 
and  visualize  the  resulting  curves. 

We  have  therefore  made  use  of  a  new  class  of  integration  algorithms  de- 
veloped by  Crouch  and  Grossman  [8],  and  Crouch,  Yan  and  Grossman  [7] 


which  do  indeed  preserve  such  structures.  The  algorithms  are  therefore  called 
geometrically  exact  [6].  We  briefly  indicate  the  important  aspects  of  these 
geometrically  stable  algorithms  which  pertain  to  the  equations  (20)  and  (21). 
Suppose  that  we  wish  to  numerically  integrate  an  ordinary  differential  equa- 
tion on  R"  given  by  the  equations: 

i(t)  =  F(t,x(t)),     a-elT,      x(0)  =  a-o  (23) 

where 

N 

F{t,x)  =  J2aJ{^x)AAx)  (24) 

n  >  A'*,  A3  are  vector  fields  on  R"  and  a°  are  functions  on  R  x  Rn.  Suppose 
that  in  addition  we  are  given  a  set  of  functions  on  Rn  whose  numerical  values 
are  constant  along  solutions  of  the  equation  (23)  and  that  the  level  sets  of 
these  functions  are  manifolds.  Denote  the  level  set  through  Xq  by  M .  It  is 
convenient  to  assume  the  slightly  stronger  assumption  that  the  vector  fields 
Aj  are  everywhere  tangent  to  M .  We  also  assume  that  there  is  an  oracle  that 
can  integrate  any  vector  field  of  the  following  form,  to  any  desired  accuracy: 

Z(x)  =  J2<*JMx)-  (25) 

i=i 

Here  aJ  are  real  numbers.  We  define  vertor  fields  Fp  by  setting: 

and  note  that  Fp  is  simply  the  vector  field  F  with  coefficients  "frozen"  at 
the  point  p.  If  we  denote  the  flow  of  any  vector  field  Z  by  (t,x)  — >  6z(t,x), 
R  x  Rn  — >  R",  then  since  the  vector  fields  Aj  are  everywhere  tangent  to  M , 
it  follows  that  x  E  M  implies  that  6pP(t,x)  €  M  for  all  p  and  for  all  t  for 
which  the  flow  is  defined. 

We  now  introduce  the  (explicit)  geometrically  exact  Runge  Kutta  algo- 
rithms as  described  in  [8].  Let  Xk  =  x(tk)  be  a  point  of  the  integral  curve 
x  of  (23).  Then,  define  vector  fields  on  Rn  by  freezing  coefficients  of  F  at 
various  points  as  follows: 

N 

Fi(x)  =  J2aJ(<t^xk)Aj(x) 

3=1 

N 

F2{x)  =  YL^^k  +  hc2i,0Fl(hc2i,xk))Aj(x) 

j=l 

8 


F3{x)  =  ]>>''(**  +  h{c31  +  cZ2),eF2{hc^OFl(hczup)))Aj(x), 

3=1 

etc.    Second,  we  define  the  numerical  integration  algorithm  via  an  update 
rule: 

z*+i  =  0Fr{hcT,9Fr_:(hcr_u---  ,6Fl{hcuxk)))  (26) 

where  h  is  the  "step  length"  and  c,  and  ctj  are  constants  to  be  determined. 
'These  constants  are  determined  from  the  "consistency  equations,"  obtained 
by  making  the  Taylor  expansions  in  /i,  about  h  =  0,  of  both  sides  of  (26), 
using  on  the  left  hand  side  the  expression  xk+i  =  9F(h^xk).  If  the  coefficients 
of  hx  agree  up  to  i  =  (/,  then  q  is  said  to  be  order  of  the  resulting  algorithm. 
Note  that  in  general  we  have  r  >  q,  while  for  classical  Runge  Kutta  schemes 
we  can  always  take  r  =  q.  Note  that  the  update  rule  defined  by  equation 
(26)  has  the  property  that  if  xk  lies  in  A/,  then  so  does  xk+i  since  each  flow  is 
defined  by  a  vector  field  F  with  frozen  coefficients.  In  the  special  case  where 
?/  —  N,  M  =  R"  and  Ax  =  e,  is  the  standard  ith  basis  vector  in  Rn ,  the 
algorithm  reduces  to  the  form  of  a  classical  explicit  Runge-Kutta  algorithm. 
In  the  paper  [8]  the  consistency  equations  are  derived  for  the  geometri- 
cally exact  Runge  Kutta  algorithms  via  a  careful  geometric  analysis  of  the 
equations  (26).  The  results  show  that  a  geometrically  exact  third  order  ex- 
plicit Runge  Kutta  algorithm  can  be  obtained  for  r  —  3  and  is  determined 
by  five  independent  consistency  equations,  in  the  six  constants  defining  the 
algorithm.  The  equations  have  multiple  solutions,  all  of  which  are  solutions 
to  the  equations  which  determine  the  classical  explicit  Runge  Kutta  algo- 
rithms. Thus,  all  solutions  define  classical  algorithms,  but  the  solutions  are 
not  ones  traditionally  found  in  the  literature.  We  have  used  the  following 
solution  of  all  five  equations: 

Cl  =  l,     c2  =  -2/3,     c3  =  2/3, 

c21  =  -l/24,      c31  =  161/24,      c32  =  -6. 

In  the  special  case  of  the  set  of  equations  (20)  and  (21),  defining  the 
elastica  on  a  sphere,  we  use  a  hybrid  algorithm  in  which  equations  (20)  are 
integrated  using  the  third  order  geometrically  exact  Runge  Kutta  algorithm 
described  above  by  freezing  the  coefficients  Kg.  These  coefficients  are  then 
updated  by  integrating  equation  (21)  using  elliptic  functions. 

In  equation  (20),  M  is  the  three  dimensional  submanifold  of  R9  deter- 
mined by  the  equations 

<  X\jXi  >=  r2,      <  Xi,x2  >=  0,      <  Xx,x3  >=  0, 

<22,22>=1,        <  22,^3  >=0,        <ar3,a;3>=l. 


The  solutions  of  (20)  for  any  initial  condition  x(Q)  =  x0  £  M  lie  completely 
in  M.  This  follows  from  the  fact  that  the  functions 

fi=<xuXi>,     f2=<xux2>,     h=<xx,xz>, 

f4   =<  X2,X2  >i        J5  —  <  x2tX3  >,        /6  =  <  X3,  X3  > 

satisfy  the  system  of  differential  equations 

f[     =    2/a 

ti     =    /4  +  s/3-(l/r2)/! 

fz  =  h  —  Kgh 

f'<  =  2n9h  -  (2/r2)/2 

/s  =  ^fe-(l/r2)f3-K9f4 

fe  =  -2«,/5 

which  has  the  unique  solution 

/i  =  r2,  /a  =  0,  f3  =  0,  /,  =  1,  fB  =  0,  /6  =  1. 

Since  this  argument  does  not  specify  the  function  ac5,  it  follows  that  the  flows 
of  the  vector  fields  F\  —  Fz  with  frozen  coefficients  are  mappings  into  M  as 
needed  in  formula  (26). 

Furthermore,  the  vector  fields  Aj  in  equation  (20)  are  linear,  so  that  F  is 
given  by  an  expression  of  the  form 

F(x)  =  '£V(tix)Bix  (27) 

3=1 

where  Bj  are  matrices  and  b1  are  functions.  Freezing  the  functions  b3  to  values 
f33  yields  a  system  of  linear  differential  equations  with  constant  coefficients: 

Thus,  the  flow  9Fi  of  the  vector  field  F,  (i= 1,2,3)  is  given  by 

^F.C?)  =exp(td)q 
and  (26)  takes  the  special  form 

xk+i  —  exp(c3/iC3)  '  ex.\>{c2hC2)  ■  exp(cihC\)  •  x^. 
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Since  61  =  \/r  is  constant  and  b2  —  Kg  depends  only  on  t,  the  matrices  Cx 
are  given  by 

Ci  =  C(tk),     C2  =  C(tk  +  hc2l),     C3  =  C(tk  +  h(c31  +  c3i))        I 

where  C(t)  :=  12j=i  V{i)Bj.  C  can  be  considered  as  a  3  x  3  skew  symmetric 
matrix  with  matrix  components  that  are  themselves  3x3  matrices.  Therefore 
the  standard  formula  for  the  exponential  of  a  3  x  3  skew  symmetric  matrix 
can  be  used  to  compute  the  flow  of  each  vector  field: 

exp{t(f)S{c))  =  I  +  sm{t<j>)S{c)  +  (1  -  cos{t<f>))S{c)2 

where  S(c)  is  the  skew  symmetric  matrix  satisfying  S(a)b  =  bxa,  and  \c\  =  1. 

Regarding  the  performance  of  the  geometrically  stable  integration  method 
it  is  clear  that  this  method  is  computationally  more  expensive  than  a  classical 
algorithm  with  the  same  number  of  stages  and  step  length.  The  exact  cost 
of  the  new  integration  scheme  can  be  found  in  [6].  However,  if  we  compare 
the  performance  of  the  geometrically  stable  method  with  the  classical  fourth 
order  Runge-Kutta  algorithm  on  the  problem  of  spherical  elastica,  it  turns 
out  that  a  slightly  increased  step  length  for  the  new  method  suffices  to  out- 
perform the  classical  integration  method.  Using  MATLAB  implementations 
of  both  algorithms  we  found  that  the  geometrically  stable  method  needs  an 
average  value  of  1001  Mflops  for  one  complete  step  while  the  fourth  order 
Runge-Kutta  only  uses  an  average  value  of  768  Mflops  per  step.  Therefore, 
one  can  statistically  achieve  the  same  performance  by  using  an  increased  step 
length  of  h  «  1.3038  /i/?a'  f°r  the  new  method.  Taking  into  account  that  the 
proposed  algorithm  not  only  delivers  points  that  lie  exactely  on  the  sphere 
but  that  are  also  approximately  equally  spaced  along  the  tracked  curve,  the 
geometrically  stable  method  seems  to  be  a  good  choice  for  the  integration  of 
spherical  elastica. 

A  performance  comparison  of  the  new  integration  scheme  with  a  more 
sophisticated  classical  method,  the  IMSL  Runge  Kutta  implementation,  can 
be  found  in  [71. 


5      Classification  of  spherical  elastica 

Acting  on  the  suggestion  of  D.  Bernoulli,  L.  Euler  derived  differential  equa- 
tions for  plane  elastica  and  classified  the  fundamental  forms  of  these  curves 
(see  [9],  [11]).  A  curvature  analysis  of  the  various  fundamental  cases  has  been 
given  in  [3]. 

In  this  section  we  classify  the  forms  of  spherical  elastica  based  on  the 
differential  equation  (21)  for  the  geodesic  curvature.  This  equation  is  of  the 
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same  form  as  the  equation  for  the  curvature  of  plane  elastica  and  can  be 
solved  in  terms  of  Jacobi's  elliptic  functions  in  the  form: 

Kg(s)  =  Km  dn(«m(s  -  sm)/2  1 12) 

where  the  positive  parameter  I2  of  the  elliptic  function  is  given  by 

t>  =  2(K-+f2-g)  (29) 

(see  [3]).  The  parameter  Km  represents  the  amplitude  of  the  periodic  curva- 
ture function  and  sm  denotes  the  value  at  which  Ac(sm)  =  Km. 

To  obtain  a  representation  of  the  curvature  in  terms  of  Jacobi's  functions 
with  parameter  /2  smaller  than  1,  one  uses  the  formula 


«,M  =  Km  cn{yJ(Kl  +  2/r*-a)/2(s  -  sm)  I  i)  (30) 

if  <j  <  0.5/c2^  +  2/r2.  Since  the  function  en  has  zeros  while  dn  is  positive,  the 
above  case  distinction  reflects  the  main  division  of  elastic  curves  into  those 
where  the  geodesic  curvature  changes  sign  and  the  other  with  constant  sign 
of  their  geodesic  curvature. 

The  change  of  the  forms  of  spherical  elastica  while  Km  is  fixed  and  a  in- 
creases is  shown  by  figures  1-8.  The  maximum  value  of  the  tension  parameter 
a  for  a  real  elastic  curve  on  a  sphere  is  according  to  (29)  a  =  Km  +  2/r2. 
This  choice  of  a  corresponds  to  the  dashed  circle  that  is  shown  in  all  figures 
for  the  purpose  of  orientation.  The  second  curve  in  figure  1  has  a  negative 
tension  parameter  of  high  absolute  value  (<r  =  —10000).  In  comparison  to 
figure  2  where  a  =  —30  we  observe  that  decreasing  the  tension  parameter 
has  the  effect  of  lowering  the  amplitude  of  the  curve  as  it  is  known  from 
the  euclidean  case  (see  [3]).  The  oscillation  of  the  curve  in  figure  1  is  of 
too  low  amplitude  to  be  visually  observed  and  the  curve  can  not  be  distin- 
guished from  the  geodesic  determined  by  the  initial  conditions.  (Note,  that 
the  geodesic  curvature  (30)  is  far  from  getting  flat  for  a  —*  — oo  but  in  fact 
approaches  a  cos  function  that  oscillates  with  a  period  that  decreases  with 
a.) 

A  curve  with  positive  tension  parameter  is  shown  in  figure  3  where  a  = 
2/r2  =  2.  The  displayed  curve  is  a  special  case  because  the  parameter 
l//2  =  1/2.  Here  the  geodesic  curvature  is  given  by  the  lemniscate  function: 

Kg(s)  =  Km  coslemn(>cm(s  -  sm)/2). 

In  figures  4  and  5  it  is  illustrated  that  with  increasing  positive  a  the 
bays  of  the  curves  start  to  overlap  until  a  figure-8-configuration  is  reached 
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Figure  1:  Km  =  4,  a  =  -10000 
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where  all  double  points  of  the  curve  coincide.  This  happens  for  a  parameter 
a  %  7.849. 

While  a  increases  further  the  curve  proceeds  through  the  figure-8-shape 
forming  a  series  of  loops  with  alternating  sign  of  geodesic  curvature  (see 
figure  6).  These  loops  recede  from  each  other  until  in  the  limiting  case,  when 
a  =  2/r2  +  0.5k^,  the  curve  forms  a  single  loop  (see  figure  7).  Here  the 
geodesic  curvature  is  given  by 

Kg{s)  -  Kmsech{Km(s  -  sm)/2). 

Figure  8  shows  that  the  single  loop  transforms  into  a  series  of  loops  with 
the  same  sign  of  geodesic  curvature.  With  increasing  a  the  loops  come  closer 
together  and  finally  collapse  into  a  circle  when  a  =  2/r    -f  k^. 
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